# Code to create Figure A10 in "Pandemics and Cities: Evidence from the Black 
# Death and the Long-Run"
# Author: Noel Johnson
# Date Created: 10-17-19
# Last updated: 9-8-23

library(raster)
library(sf)
library(tidyverse)
library(tmap)

setwd("/Users/noeljohnson_laptop/Dropbox/Research/Plague City Growth/JUE RR1 Nov 2022/figures_replication/Figure_A10")

# Load in the data sets, make spatial, and project
# Define the projection
EEC <- "+proj=eqdc +lat_0=0 +lon_0=0 +lat_1=43 +lat_2=62 +x_0=0 +y_0=0 +ellps=intl +units=m +no_defs"

# Load in the cities data
cities_full <- read.csv("list_467_165_cities.csv")
head(cities_full)
# spatialize the cities data
cities_full_spatial <- st_as_sf(cities_full,
                                coords = c("longitude", "latitude"), crs = 4326)
# project the cities data
cities_165_proj <- st_transform(cities_full_spatial, EEC)
st_geometry(cities_165_proj)

# Load in the modern country borders
modern_borders <- st_read("ModernCountries/Modern Europe Projected.shp")
# project the modern country borders
modern_borders_proj <- st_transform(modern_borders, EEC)
st_geometry(modern_borders_proj)

# need to crop the maps
# original Extent of modern countries
# c(xmin=-2639987, xmax=1655719, ymin=3243609, ymax=8097861))
modern_borders_clip <- st_crop(modern_borders_proj,
                               c(xmin=-1100000, xmax=1655719,
                                 ymin=3243609, ymax=6800000))

# Make a simple plot of the cities and borders to make sure it looks right
plot(modern_borders_clip["name"], col=NA, reset = FALSE)
plot(cities_165_proj[5], add = TRUE)

# Create two data sets:
# (1) 467 cities with pop >1000 in 1300;
# (2) 165 cities in our sample (mort = 1)
head(cities_165_proj)
pop1300_cities <- cities_165_proj %>%
  select(-samplemort)
sample_cities <- cities_165_proj %>%
  filter(samplemort==1) %>%
  select(-samplemort)

# Make the figure for the paper
state_map <- tm_shape(modern_borders_clip) + tm_borders(col = "black")

breaks <- c(50, 100, 150, 200, 250)
cities_map <- tm_shape(pop1300_cities) +
  tm_bubbles("pop1300", "black", border.col = "white",
             shape = 1,
             sizes.legend = breaks,  
             title.size ="Population in 1300",
             scale = 1.75, perceptual = F, sizes.legend.labels = c("50", "100", "150", "200", "250")) +
  tm_layout(legend.stack = "vertical", legend.text.size = .75, frame = F, inner.margins = .10, legend.position = c(-.05, 0.36)) +
  state_map

full_map_x <- tm_shape(sample_cities) +
  tm_dots("black", size = .4, shape = 4) +
  tm_add_legend(col = "black", type = "symbol", shape = 4, title = "Sample Cities") +
  cities_map

full_map_x

# save the map
tmap_save(full_map_x, filename="figure_a10.png",
          height=8.5, width=11, units="in", dpi=300)


#











# End Code